/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package bd.org.apache.commons.math.analysis.solvers;

import bd.org.apache.commons.math.analysis.UnivariateFunction;
import bd.org.apache.commons.math.exception.NoBracketingException;
import bd.org.apache.commons.math.exception.NotStrictlyPositiveException;
import bd.org.apache.commons.math.exception.NullArgumentException;
import bd.org.apache.commons.math.exception.NumberIsTooLargeException;
import bd.org.apache.commons.math.exception.util.LocalizedFormats;
import bd.org.apache.commons.math.util.FastMath;

/**
 * Utility routines for {@link UnivariateSolver} objects.
 *
 * @version $Id: UnivariateSolverUtils.java 1244107 2012-02-14 16:17:55Z erans $
 */
public class UnivariateSolverUtils {

    /**
     * Class contains only static methods.
     */
    private UnivariateSolverUtils() {}

    /**
     * Convenience method to find a zero of a univariate real function.  A default
     * solver is used.
     *
     * @param function Function.
     * @param x0 Lower bound for the interval.
     * @param x1 Upper bound for the interval.
     * @return a value where the function is zero.
     */
    public static double solve(UnivariateFunction function, double x0, double x1) {

        if (function == null) {
            throw new NullArgumentException(LocalizedFormats.FUNCTION);
        }

        final UnivariateSolver solver = new BrentSolver();

        return solver.solve(Integer.MAX_VALUE, function, x0, x1);
    }


    /**
     * Convenience method to find a zero of a univariate real function.  A default
     * solver is used.
     *
     * @param function Function.
     * @param x0 Lower bound for the interval.
     * @param x1 Upper bound for the interval.
     * @param absoluteAccuracy Accuracy to be used by the solver.
     * @return a value where the function is zero.
     */
    public static double solve(UnivariateFunction function, double x0, double x1,
                               double absoluteAccuracy) {

        if (function == null) {
            throw new NullArgumentException(LocalizedFormats.FUNCTION);
        }

        final UnivariateSolver solver = new BrentSolver(absoluteAccuracy);

        return solver.solve(Integer.MAX_VALUE, function, x0, x1);
    }


    /**
     * Force a root found by a non-bracketing solver to lie on a specified side,
     * as if the solver was a bracketing one.
     *
     * @param maxEval maximal number of new evaluations of the function
     * (evaluations already done for finding the root should have already been subtracted
     * from this number)
     * @param f function to solve
     * @param bracketing bracketing solver to use for shifting the root
     * @param baseRoot original root found by a previous non-bracketing solver
     * @param min minimal bound of the search interval
     * @param max maximal bound of the search interval
     * @param allowedSolution the kind of solutions that the root-finding algorithm may
     * accept as solutions.
     * @return a root approximation, on the specified side of the exact root
     */
    public static double forceSide(final int maxEval, final UnivariateFunction f,
                                   final BracketedUnivariateSolver<UnivariateFunction> bracketing,
                                   final double baseRoot, final double min, final double max,
                                   final AllowedSolution allowedSolution) {

        if (allowedSolution == AllowedSolution.ANY_SIDE) {

            // no further bracketing required
            return baseRoot;
        }

        // find a very small interval bracketing the root
        final double step = FastMath.max(bracketing.getAbsoluteAccuracy(),
                                         FastMath.abs(baseRoot * bracketing.getRelativeAccuracy()));
        double xLo           = FastMath.max(min, baseRoot - step);
        double fLo           = f.value(xLo);
        double xHi           = FastMath.min(max, baseRoot + step);
        double fHi           = f.value(xHi);
        int    remainingEval = maxEval - 2;

        while (remainingEval > 0) {

            if (((fLo >= 0) && (fHi <= 0)) || ((fLo <= 0) && (fHi >= 0))) {

                // compute the root on the selected side
                return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution);
            }

            // try increasing the interval
            boolean changeLo = false;
            boolean changeHi = false;

            if (fLo < fHi) {

                // increasing function
                if (fLo >= 0) {
                    changeLo = true;
                } else {
                    changeHi = true;
                }
            } else if (fLo > fHi) {

                // decreasing function
                if (fLo <= 0) {
                    changeLo = true;
                } else {
                    changeHi = true;
                }
            } else {

                // unknown variation
                changeLo = true;
                changeHi = true;
            }

            // update the lower bound
            if (changeLo) {

                xLo = FastMath.max(min, xLo - step);
                fLo = f.value(xLo);

                remainingEval--;
            }

            // update the higher bound
            if (changeHi) {

                xHi = FastMath.min(max, xHi + step);
                fHi = f.value(xHi);

                remainingEval--;
            }

        }

        throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING, xLo, xHi, fLo, fHi,
                                        maxEval - remainingEval, maxEval, baseRoot, min, max);

    }


    /**
     * This method attempts to find two values a and b satisfying <ul>
     * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
     * <li> <code> f(a) * f(b) < 0 </code></li>
     * </ul>
     * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
     * and <code>b</code> bracket a root of f.
     * <p>
     * The algorithm starts by setting
     * <code>a := initial -1; b := initial +1,</code> examines the value of the
     * function at <code>a</code> and <code>b</code> and keeps moving
     * the endpoints out by one unit each time through a loop that terminates
     * when one of the following happens: <ul>
     * <li> <code> f(a) * f(b) < 0 </code> --  success!</li>
     * <li> <code> a = lower </code> and <code> b = upper</code>
     * -- NoBracketingException </li>
     * <li> <code> Integer.MAX_VALUE</code> iterations elapse
     * -- NoBracketingException </li>
     * </ul></p>
     * <p>
     * <strong>Note: </strong> this method can take
     * <code>Integer.MAX_VALUE</code> iterations to throw a
     * <code>ConvergenceException.</code>  Unless you are confident that there
     * is a root between <code>lowerBound</code> and <code>upperBound</code>
     * near <code>initial,</code> it is better to use
     * {@link #bracket(UnivariateFunction, double, double, double, int)},
     * explicitly specifying the maximum number of iterations.</p>
     *
     * @param function Function.
     * @param initial Initial midpoint of interval being expanded to
     * bracket a root.
     * @param lowerBound Lower bound (a is never lower than this value)
     * @param upperBound Upper bound (b never is greater than this
     * value).
     * @return a two-element array holding a and b.
     */
    public static double[] bracket(UnivariateFunction function, double initial, double lowerBound,
                                   double upperBound) {

        return bracket(function, initial, lowerBound, upperBound, Integer.MAX_VALUE);
    }


    /**
     * This method attempts to find two values a and b satisfying <ul>
     * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
     * <li> <code> f(a) * f(b) <= 0 </code> </li>
     * </ul>
     * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
     * and <code>b</code> bracket a root of f.
     * <p>
     * The algorithm starts by setting
     * <code>a := initial -1; b := initial +1,</code> examines the value of the
     * function at <code>a</code> and <code>b</code> and keeps moving
     * the endpoints out by one unit each time through a loop that terminates
     * when one of the following happens: <ul>
     * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
     * <li> <code> a = lower </code> and <code> b = upper</code>
     * -- NoBracketingException </li>
     * <li> <code> maximumIterations</code> iterations elapse
     * -- NoBracketingException </li></ul></p>
     *
     * @param function Function.
     * @param initial Initial midpoint of interval being expanded to
     * bracket a root.
     * @param lowerBound Lower bound (a is never lower than this value).
     * @param upperBound Upper bound (b never is greater than this
     * value).
     * @param maximumIterations Maximum number of iterations to perform
     * @return a two element array holding a and b.
     */
    public static double[] bracket(UnivariateFunction function, double initial, double lowerBound,
                                   double upperBound, int maximumIterations) {

        if (function == null) {
            throw new NullArgumentException(LocalizedFormats.FUNCTION);
        }

        if (maximumIterations <= 0) {

            throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS,
                    maximumIterations);
        }

        verifySequence(lowerBound, initial, upperBound);

        double a = initial;
        double b = initial;
        double fa;
        double fb;
        int    numIterations = 0;

        do {

            a  = FastMath.max(a - 1.0, lowerBound);
            b  = FastMath.min(b + 1.0, upperBound);
            fa = function.value(a);

            fb = function.value(b);

            ++numIterations;

        } while ((fa * fb > 0.0) && (numIterations < maximumIterations)
                 && ((a > lowerBound) || (b < upperBound)));

        if (fa * fb > 0.0) {

            throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING, a, b, fa, fb,
                                            numIterations, maximumIterations, initial, lowerBound,
                                            upperBound);
        }

        return new double[] { a, b };
    }


    /**
     * Compute the midpoint of two values.
     *
     * @param a first value.
     * @param b second value.
     * @return the midpoint.
     */
    public static double midpoint(double a, double b) {

        return (a + b) * 0.5;
    }


    /**
     * Check whether the interval bounds bracket a root. That is, if the
     * values at the endpoints are not equal to zero, then the function takes
     * opposite signs at the endpoints.
     *
     * @param function Function.
     * @param lower Lower endpoint.
     * @param upper Upper endpoint.
     * @return {@code true} if the function values have opposite signs at the
     *         given points.
     */
    public static boolean isBracketing(UnivariateFunction function, final double lower,
                                       final double upper) {

        if (function == null) {
            throw new NullArgumentException(LocalizedFormats.FUNCTION);
        }

        final double fLo = function.value(lower);
        final double fHi = function.value(upper);

        return ((fLo >= 0) && (fHi <= 0)) || ((fLo <= 0) && (fHi >= 0));
    }


    /**
     * Check whether the arguments form a (strictly) increasing sequence.
     *
     * @param start First number.
     * @param mid Second number.
     * @param end Third number.
     * @return {@code true} if the arguments form an increasing sequence.
     */
    public static boolean isSequence(final double start, final double mid, final double end) {

        return (start < mid) && (mid < end);
    }


    /**
     * Check that the endpoints specify an interval.
     *
     * @param lower Lower endpoint.
     * @param upper Upper endpoint.
     */
    public static void verifyInterval(final double lower, final double upper) {

        if (lower >= upper) {

            throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL, lower,
                    upper, false);
        }
    }


    /**
     * Check that {@code lower < initial < upper}.
     *
     * @param lower Lower endpoint.
     * @param initial Initial value.
     * @param upper Upper endpoint.
     * {@code initial >= upper}.
     */
    public static void verifySequence(final double lower, final double initial,
                                      final double upper) {

        verifyInterval(lower, initial);
        verifyInterval(initial, upper);
    }


    /**
     * Check that the endpoints specify an interval and the end points
     * bracket a root.
     *
     * @param function Function.
     * @param lower Lower endpoint.
     * @param upper Upper endpoint.
     */
    public static void verifyBracketing(UnivariateFunction function, final double lower,
            final double upper) {

        if (function == null) {
            throw new NullArgumentException(LocalizedFormats.FUNCTION);
        }

        verifyInterval(lower, upper);

        if (!isBracketing(function, lower, upper)) {

            throw new NoBracketingException(lower, upper, function.value(lower),
                                            function.value(upper));
        }
    }
}
